/**
 * @file  hydrocode.c
 * @brief This is a C file of the main function.
 */

/** 
 * @mainpage 2D Godunov/GRP scheme for Eulerian hydrodynamics
 * @brief This is an implementation of fully explict forward Euler scheme
 *        for 2-D Euler equations of motion on Eulerian coordinate.
 * @version 0.2
 *
 * @section File_directories File directories
 * <table>
 * <tr><th> data_in/  <td> Folder to store input files RHO/U/P/config.txt
 * <tr><th> data_out/ <td> Folder to store output files RHO/U/P/E/X/log.txt
 * <tr><th> doc/      <td> Code documentation generated by doxygen
 * <tr><th> src/      <td> Folder to store C source code
 * </table>
 * 
 * @section Program_structure Program structure
 * <table>
 * <tr><th> include/                   <td> Header files
 * <tr><th> tools/                     <td> Tool functions
 * <tr><th> file_io/                   <td> Program reads and writes files
 * <tr><th> riemann_solver/            <td> Riemann solver programs
 * <tr><th> inter_process/             <td> Intermediate processes in finite volume scheme
 * <tr><th> flux_calc/                 <td> Program for calculating numerical fluxes in finite volume scheme
 * <tr><th> finite_volume/             <td> Finite volume scheme programs
 * <tr><th> hydrocode_2D/hydrocode.c   <td> Main program
 * <tr><th> hydrocode_2D/hydrocode.sh  <td> Bash script compiles and runs programs
 * </table>
 *
 * @section Exit_status Program exit status code
 * <table>
 * <tr><th> exit(0)  <td> EXIT_SUCCESS
 * <tr><th> exit(1)  <td> File directory error
 * <tr><th> exit(2)  <td> Data reading/writing error
 * <tr><th> exit(3)  <td> Calculation error
 * <tr><th> exit(4)  <td> Arguments error
 * <tr><th> exit(5)  <td> Memory error
 * </table>
 * 
 * @section Compile_environment Compile environment
 *          - Linux/Unix: gcc, glibc, MATLAB/Octave
 *            - Compile in 'src/hydrocode_2D': Run './hydrocode.sh' command on the terminal.
 *          - Winodws: Visual Studio, MATLAB/Octave
 *            - Create a C++ Project from Existing Code in 'src/hydrocode_2D/' with ProjectName 'hydrocode'.
 *            - Compile in 'x64/Debug' using shortcut key 'Ctrl+B' with Visual Studio.
 *
 * @section Usage_description Usage description
 *          - Input files are stored in folder 'data_in/two-dim/name_of_test_example/'.
 *          - Input files may be produced by MATLAB/Octave script 'value_start.m'.
 *          - Description of configuration file 'config.txt/.dat' refers to 'doc/config.csv'.
 *          - Run program:
 *            - Linux/Unix: Run 'shell/hydrocode_run.sh' command on the terminal. \n
 *                          The details are as follows: \n
 *                          Run 'hydrocode.out name_of_test_example name_of_numeric_result order[_scheme]
 *                               coordinate config[n]=(double)C' command on the terminal. \n
 *                          e.g. 'hydrocode.out GRP_Book/6_1 GRP_Book/6_1 2[_GRP] EUL 5=100' (second-order Eulerian GRP scheme).
 *                          - order: Order of numerical scheme (= 1 or 2).
 *                          - scheme: Scheme name (= Riemann_exact/Godunov, GRP or …).
 *                          - coordinate: Eulerian coordinate framework (= EUL).
 *            - Windows: Run 'hydrocode.bat' command on the terminal. \n
 *                       The details are as follows: \n
 *                       Run 'hydrocode.exe name_of_test_example name_of_numeric_result order[_scheme] 
 *                            coordinate n=C' command on the terminal. \n
 *                       [Debug] Project -> Properties -> Configuration Properties -> Debugging \n
 *             <table>
 *             <tr><th> Command Arguments <td> name_of_test_example name_of_numeric_result order[_scheme] coordinate n=C
 *             <tr><th> Working Directory <td> hydrocode_2D
 *             </table>
 *                       [Run] Project -> Properties -> Configuration Properties -> Linker -> System \n
 *             <table>
 *             <tr><th> Subsystem <td> (/SUBSYSTEM:CONSOLE)
 *             </table>
 *                       [Run] Project -> Properties -> Configuration Properties -> C/C++ -> Language \n
 *             <table>
 *             <tr><th> OpenMP Support <td> (/openmp)
 *             </table>
 * 
 *          - Output files can be found in folder 'data_out/two-dim/'.
 *          - Output files may be visualized by MATLAB/Octave script 'value_plot.m'.
 * 
 * @section Precompiler_options Precompiler options
 *          - NODATPLOT: in hydrocode.c. (Default: undef)
 *          - NOTECPLOT: in hydrocode.c. (Default: undef)
 *          - HDF5PLOT:  in hydrocode.c. (Default: undef)
 *          - EXACT_TANGENT_DERIVATIVE: in linear_GRP_solver_Edir_G2D.c.   (Default: undef)
 *          - Riemann_solver_exact_single: in riemann_solver.h.            (Default: Riemann_solver_exact_Ben)
 *          - MULTIFLUID_BASICS: 'Switch whether to compute multi-fluids.' (Default: undef)
 */


#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "../include/var_struc.h"
#include "../include/file_io.h"
#include "../include/finite_volume.h"


#ifdef DOXYGEN_PREDEFINED
/**
 * @def NODATPLOT
 * @brief Switch whether to plot without Matrix data.
 */
#define NODATPLOT
/**
 * @def HDF5PLOT
 * @brief Switch whether to plot with HDF5 data.
 */
#define HDF5PLOT
/**
 * @def NOTECPLOT
 * @brief Switch whether to plot without Tecplot data.
 */
#define NOTECPLOT
#endif

double config[N_CONF]; //!< Initial configuration data array.

/**
 * @brief N memory allocations to the initial fluid variable 'v' in the structure cell_var_stru.
 */
#define CV_INIT_MEM(v, N)						\
    do {								\
    for(k = 0; k < N; ++k)						\
	{								\
	    CV[k].v = (double **)malloc(n_x * sizeof(double *));	\
	    if(CV[k].v == NULL)						\
		{							\
		    printf("NOT enough memory! CV[%d].%s\n", k, #v);	\
		    retval = 5;						\
		    goto return_NULL;					\
		}							\
	    for(j = 0; j < n_x; ++j)					\
		{							\
		    CV[k].v[j] = (double *)malloc(n_y * sizeof(double)); \
		    if(CV[k].v[j] == NULL)				\
			{						\
			    printf("NOT enough memory! CV[%d].%s[%d]\n", k, #v, j); \
			    retval = 5;					\
			    goto return_NULL;				\
			}						\
		}							\
	}								\
    } while (0)

/**
 * @brief This is the main function which constructs the
 *        main structure of the 2-D Eulerian hydrocode.
 * @param[in] argc: ARGument Counter.
 * @param[in] argv: ARGument Values.
 *            - argv[1]: Folder name of test example (input path).
 *            - argv[2]: Folder name of numerical results (output path).
 *            - argv[3]: Order of numerical scheme[_scheme name] (= 1[_Riemann_exact] or 2[_GRP]).
 *            - argv[4]: Eulerian coordinate framework (= EUL).
 *            - argv[5,6,…]: Configuration supplement config[n]=(double)C (= n=C).
 * @return Program exit status code.
 */
int main(int argc, char *argv[])
{
  int k, i, j, retval = 0;
  // Initialize configuration data array
  for(k = 1; k < N_CONF; k++)
      config[k] = INFINITY;
  char * scheme = NULL; // Riemann_exact(Godunov), GRP
  arg_preprocess(4, argc, argv, scheme);

  // Set dimension.
  config[0] = (double)2; // Dimensionality = 2

  // The number of times steps of the fluid data stored for plotting.
  int N, N_plot;
  double * time_plot;
    /* 
     * We read the initial data files.
     * The function initialize return a point pointing to the position
     * of a block of memory consisting (n_x*n_y) variables of type double.
     * The (n_x*n_y) array elements of these variables are the initial value.
     */
  struct flu_var FV0 = initialize_2D(argv[1], &N, &N_plot, &time_plot); // Structure of initial data array pointer.
    /* 
     * (n_x*n_y) is the number of initial value as well as the number of grids.
     * As (n_x*n_y) is frequently use to represent the number of grids,
     * we do not use the name such as num_cell here to correspond to
     * notation in the math theory.
     */
  const int n_x = (int)config[13], n_y = (int)config[14];
  const double h_x = config[10], h_y = config[11], gamma = config[6];
  const int order = (int)config[9];
  const _Bool dim_split = (_Bool)config[33]; // Dimensional splitting?

  // Structure of fluid variables in computational cells array pointer.
  struct cell_var_stru * CV = (struct cell_var_stru*)malloc(N * sizeof(struct cell_var_stru));
  double ** X, ** Y;
  double * cpu_time = (double *)malloc(N * sizeof(double));
  X = (double **)malloc((n_x+1) * sizeof(double *));
  Y = (double **)malloc((n_x+1) * sizeof(double *));
  if(cpu_time == NULL)
      {
	  printf("NOT enough memory! CPU_time\n");
	  retval = 5;
	  goto return_NULL;
      }
  if(X == NULL || Y == NULL)
      {
	  printf("NOT enough memory! X or Y\n");
	  retval = 5;
	  goto return_NULL;
      }
  for(j = 0; j <= n_x; ++j)
  {
    X[j] = (double *)malloc((n_y+1) * sizeof(double));
    Y[j] = (double *)malloc((n_y+1) * sizeof(double));
    if(X[j] == NULL || Y[j] == NULL)
    {
      printf("NOT enough memory! X[%d] or Y[%d]\n", j, j);
      retval = 5;
      goto return_NULL;
    }
  }
  if(CV == NULL)
      {
	  printf("NOT enough memory! Cell Variables\n");
	  retval = 5;
	  goto return_NULL;
      }
  // Initialize arrays of fluid variables in cells.
  CV_INIT_MEM(RHO, N);
  CV_INIT_MEM(U, N);
  CV_INIT_MEM(V, N);
  CV_INIT_MEM(P, N);
  CV_INIT_MEM(E, N);
  // Initialize the values of energy in computational cells and (x,y)-coordinate of the cell interfaces.
  for(j = 0; j <= n_x; ++j)
      for(i = 0; i <= n_y; ++i)	
	  {
	      X[j][i] = j * h_x;
	      Y[j][i] = i * h_y;
	  }
  for(j = 0; j < n_x; ++j)
      for(i = 0; i < n_y; ++i)	
	  {
	      CV[0].RHO[j][i] = FV0.RHO[i*n_x + j];
	      CV[0].U[j][i]   = FV0.U[i*n_x + j];
	      CV[0].V[j][i]   = FV0.V[i*n_x + j];
	      CV[0].P[j][i]   = FV0.P[i*n_x + j];
	      CV[0].E[j][i]   = 0.5*CV[0].U[j][i]*CV[0].U[j][i] + CV[0].P[j][i]/(gamma - 1.0)/CV[0].RHO[j][i];
	      CV[0].E[j][i]  += 0.5*CV[0].V[j][i]*CV[0].V[j][i];
	  }

  if (strcmp(argv[4],"EUL") == 0) // Use GRP/Godunov scheme to solve it on Eulerian coordinate.
      {
	  config[8] = (double)0;
	  switch(order)
	      {
	      case 1:
		  config[41] = 0.0; // alpha = 0.0
	      case 2:
		  if (dim_split)
		      GRP_solver_2D_split_EUL_source(n_x, n_y, CV, X, Y, cpu_time, argv[2], N, &N_plot, time_plot);
		  else
		      GRP_solver_2D_EUL_source(n_x, n_y, CV, X, Y, cpu_time, argv[2], N, &N_plot, time_plot);
		  break;
	      default:
		  printf("NOT appropriate order of the scheme! The order is %d.\n", order);
		  retval = 4;
		  goto return_NULL;
	      }
      }
  else
      {
	  printf("NOT appropriate coordinate framework! The framework is %s.\n", argv[4]);
	  retval = 4;
	  goto return_NULL;
      }

  // Write the final data down.
#ifndef NODATPLOT
  file_2D_write(n_x, n_y, N_plot, CV, X, Y, cpu_time, argv[2], time_plot);
#endif
#ifdef HDF5PLOT
  file_2D_write_HDF5(n_x, n_y, N_plot, CV, X, Y, cpu_time, argv[2], time_plot);
#endif
#ifndef NOTECPLOT
  file_2D_write_POINT_TEC(n_x, n_y, 1, CV + N_plot-1, X, Y, cpu_time, argv[2], time_plot + N_plot-1);
#endif

 return_NULL:
  free(FV0.RHO);
  free(FV0.U);
  free(FV0.V);
  free(FV0.P);
  FV0.RHO = NULL;
  FV0.U   = NULL;
  FV0.V   = NULL;
  FV0.P   = NULL;
  for(k = 0; k < N; ++k)
  {
    for(j = 0; j < n_x; ++j)
	{
            free(CV[k].RHO[j]);
            free(CV[k].U[j]);
            free(CV[k].V[j]);
            free(CV[k].P[j]);
            free(CV[k].E[j]);
            CV[k].RHO[j] = NULL;
            CV[k].U[j]   = NULL;
            CV[k].V[j]   = NULL;
            CV[k].P[j]   = NULL;
            CV[k].E[j]   = NULL;
	}
    free(CV[k].RHO);
    free(CV[k].U);
    free(CV[k].V);
    free(CV[k].P);
    free(CV[k].E);
    CV[k].RHO = NULL;
    CV[k].U   = NULL;
    CV[k].V   = NULL;
    CV[k].P   = NULL;
    CV[k].E   = NULL;
  }
  free(CV);
  CV = NULL;
  for(j = 0; j <= n_x; ++j)
  {
      free(X[j]);
      free(Y[j]);
      X[j] = NULL;
      Y[j] = NULL;
  }
  free(X);
  free(Y);
  X = NULL;
  Y = NULL; 
  free(cpu_time);
  cpu_time = NULL;
  
  return retval;
}
